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ABSTRACT 

We study rapidly accreting, gravitationally unstable disks with a series of global, three dimensional, 
numerical experiments using the code ORION. In this paper we conduct a numerical parameter study 
focused on protostellar disks, and show that one can predict disk behavior and the multiplicity of the 
accreting star system as a function of two dimensionless parameters which compare the disk's accretion 
rate to its sound speed and orbital period. Although gravitational instabilities become strong, we find 
that fragmentation into binary or multiple systems occurs only when material falls in several times 
more rapidly than the canonical isothermal limit. The disk-to-star accretion rate is proportional to the 
infall rate, and governed by gravitational torques generated by low-m spiral modes. We also confirm 
the existence of a maximum stable disk mass: disks that exceed ~ 50% of the total system mass are 
subject to fragmentation and the subsequent formation of binary companions. 
Subject headings: star formation, accretion disks, hydrodynamics, binary stars 



1. INTRODUCTION 

Accretion disks are ubiquitous in astrophysical sys- 
tems, feeding objects ranging in size from planets to black 
holes. Rapid disk accretion requires rapid angular mo- 
mentum transport, but only a few transport mechanisms 
are known: self-gravity, turbulence from the magnetoro- 
tational instability, driven hydrodynamical turbulence, 
and (in rare circumstances) magneto-centrifugal winds. 
Our focus is on the dynamics of disks around young, 
rapidly accreti ngs protostars, for which self-gravity is the 
key ingred ient ( |Lin fc Pringle|1987 Gammie|2001 Krat 
ter et al. 2008). Gravitational instability (hereafter CI 



ility (hereafter CI) 
is important whenever disks are cold enough or massive 
enough to trigger it, as they typically are during the early 
phases of star formation. GI plays a strong role in AGN 
disks as well, and possibly in other contexts where disks 
are cold and accretion is fast. 

The role of GI in angular momentum transport is com- 
plicated by the fact that it can lead to runaway col- 
lapse. Indeed, this makes GI an attractive mechanism 
for the formation of stars with companions (binaries, 
brown dwarfs, even planets in some circumstances: |Bon- 



During the earliest phases of star formation, protostel- 
lar disks are deeply embedded within their natal clouds. 
Observing this stage has been difficult because the source 
is only visible at wavelengths where resolution is poor. 
But because it sets the initial conditions for stellar and 
planetary systems, an understanding of this phase is crit- 
ical. Models of rapidly accreting disks, like those pre- 
sented here, will be useful for the interpretation of obser- 
vations of future facilities such as ALMA and the EVLA. 

Although semi-analytical and low-dimensional studies 
can illuminate trends and provide useful approximate re- 
sults, disk fragmentation is inherently a nonlinear and 
multidimensional process. For this reason we have em- 
barked on a survey of global, three-dimensional, numer- 
ical experiments to examine the role of GI as the media- 
tor of the accretion rate in self-gravitating disks, and as 
a mechanism for creating disk-born companions. 

Any such project faces a central difficulty: the GI is 
famously sensitive to the disk's thermodynamics (Boss 



nell fc Bat e 1994a bj [Whitworth et al.||2006|). |Kratter fc 



5| |Mejia et al.|20b5||Krumholz|2006|jCai et al [ 2008 ) . 

lie it is possible and valuable to incorporate detailed 



et al.|2000[ |Gammie|200"Tl |Pickett et al 12003 1 |Rice et al. 



Matzn er 2006 (2006; hereafter KM06) and Kratte reTal 
2008| (2008,~nereafter KMK08) found that disk fragmen- 
tation and binary formation are increasingly likely as one 
considers more and more massive stars, whereas disks in 



low-mass st ar formation arc relatively stable ( Matzner & 
Levin|2005 ) . The increased frequency of giant planet for 



mation around A stars relative to F and G stars, and its 



lower sensitivity to the stellar metallicity ( Johnson et al. 
20071, may also be related. 



2005 

While it is possible and valuable to incorporate 
heating and cooling into numerical simulations as has 
been explored by the authors above, there is a cost: sim- 
ulations with these important physical processes cannot 
be scaled to represent a whole range of physical environ- 
ments, whereas those without them can. We choose to 
separate the dynamical problem from the thermal one. 
We exclude thermal physics from our simulations en- 
tirely, while scanning a thermal parameter in our survey. 
By this means we reduce the physical problem to two di- 
mensionless parameters: one for the disk's temperature, 
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another for its rotation period - both in units determined 
by its mass accretion rate. We hold these fixed in each 
simulation by choosing well-controlled initial conditions 
corresponding to self-similar core collapse. This param- 
eterization is a central aspect of our work: it forms the 
basis for our numerical survey; it allows us to treat as- 
trophysically relevant disks, including fragmentation and 
the formation of binary companions, while also maintain- 
ing generality; and it distinguishes our work from previ- 
ous numerical studies of core collapse, disk formation, 
GI, and fragmentation. 

This paper, the first in a series, focuses on the broad 
conclusions we can draw from our parameter space study; 
subsequent papers will discuss the detailed behavior of 
multiple systems, three dimensional effects such as tur- 
bulence, and vertical flows, and non-linear GI mode cou- 
pling. We begin here by introducing our dimensionless 
parameters in S|2j We describe the initial conditions and 
the numerical code used in £j3] In Q we derive analytic 
predictions for the behavior of disks as a function of our 
parameters. We describe the main results from our nu- 
merical experiments in ^5] with more detailed analysis in 
S|6] We compare them in detail to other numerical and 
analytic models of star formation in [|7] 

2. A NEW PARAMETER SPACE FOR ACCRETION 

We consider the gravitational collapse of a rotating, 
quasi-spherical gas core onto a central pointlike object, 
mediated by a disk. In the idealized picture we will ex- 
plore in this paper, the disk and the mass flows into and 
out of it can be characterized by a few simple parameters. 
At any given time, the central point mass (or masses, in 
cases where fragmentation occurs) has mass M„ , the disk 
has mass M^, and the combined mass of the two is M^. 
The disk is characterized by a constant sound speed c St d- 
Material frorn the core falls onto the disk with a mass ac- 
cretion rate M- m , and this material carries mean specific 
angular momentum (j)in, an d as a result it circularizes 
and goes into Keplerian rotation at some radius Rk,in', the 
angular velocity of the orbit is f2fc,i n - In general in what 
follows, we refer to quantities associated with the central 
object with subscript *, quantities associated with the 
disk with a subscript d, quantities associated with infall 
with subscript in. Angle brackets indicate mass-weighted 
averages over the disk (with subscript d) or over infalling 
mass (with subscript in). 

Our simple decom position of the problem is moti- 
vated by the work of |Gamrfiie (2001 1, |Matzner fc Levin 
( |2005[ ), and KMK08. We characterize our numerical ex- 
periments using two dimensionless parameters which are 
well-adapted to systems undergoing rapid accretion. We 
encapsulate the complicated physics of heating and cool- 
ing through the thermal parameter 



19931. 

The second, rotational parameter 



M in G 

„3 ' 



(1) 



which relates the infall mass accretion rate M- ln to the 
characteristic sound speed c s .d of the disk material. Our 
parameter £ is also related to the physics of core collapse 
leading to star formation. If the initial core is character- 



ized by a signal speed c e ff, c then 



JG, imply- 



ing £ ~ c off c / c sd ~ although there can be large varia- 
tions around this value ( Larson|1972 Foster & Chevalier 



r = 



M in (j)fr 



(2) 



compares the system growth rate or accretion timescale, 
M- m /M*d to the orbital timscale of infalling gas. Un- 
like £, T is independent of disk heating and cooling, de- 
pending only on the core structure and velocity field. 
In general, T compares the relative strength of rotation 
and gravity in the core. Systems with a large value of T 
(e.g. accretion-induced collapse of a white dwarf) gain 
a significant amount of mass in each orbit, and tend to 
be surrounded by thick, massive accretion disks, while 
those with very low T (e.g. active galactic nuclei) grow 
over many disk lifetimes, and tend to harbor thin disks 
with little mass relative to the central object. We con- 
sider characteristic values for our parameters in §2.1| and 
their evolution in the isothermal collapse of a rigidly ro- 
tating Bonnor-Ebert sphere in |7.1| 

In addition to being physically motivated, our param- 
eters are practical from a numerical and observational 
perspective. The thermal parameter £ is straightforward 
to calculate in other simulations, theoretical models, and 
observed disks. Accretion rates are routi nely estimated 
from measures of infall velocities in cores (Cesaroni et al. 
|2007| although some uncertainties persist}^ Disk tem- 
peratures can also be measured using infrared and sub- 
millimeter disk detections. By contrast, measuring clas- 
sic dimensionl ess disk para meters such as Toomre's Q = 
c s O/(ttGE) flToomre||l964l l can be difficult. In observed 
disks, estimating Q is challenging due to the current 
resolution and sensitivity of even the best instruments. 
While constraining disk temperatures is possible, mea- 
suring a ccurate (within a fa ctor of ~ 3) surface densities 
are not ( |Eisner et al. 2008). The rotation parameter T 
is similarly practical: one need not know the density dis- 
tribution of the initial core, nor the radial and angular 
distribution of the velocity profile - a mean value for j 
and an estimate of the core mass is sufficient to make an 
estimate for the disk size, and thus T. 

To study the evolution of systems as they accrete, we 
hold £ and T fixed for each experiment via the self -similar 
collapse of a rotating, isothermal sphere (53.2 1. This 



strategy allows us to map directly between the input pa- 
rameters, and relevant properties of the system. Specifi- 
cally, we expect dimensionless properties like the disk-to- 
star mass ratio, Toomre parameter, stellar multiplicity, 
etc., to fluctuate around well-defined mean values (see 



f3.4| 



We aim to use our parameters £ and T to: (a) explore 
the parameter space relevant to a range of star forma- 
tion scenarios; (b) better understand the disk parame- 
ters, both locally and globally, which dictate the disk 
accretion rate and fragmentation properties; (c) make 
predictions for small scale disk behavior based on larger 
scale, observable quantities; and (d) allow the results of 
more complicated and computationally expensive simu- 
lations to be extended into other regimes. 

2.1. Characteristic values of the accretion parameters 

We base our estimates of T and £ on observations 
of core rotation in low-mass and massive star-forming 
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regions (|Myers fc Fuller 1992 [Goodman et al. 

T 



Williams & Myers 1999J7 as well as the analyt ica 



1993 



esti- 



mates of cbrc rotat ion and disk te mperature in |Matzner| 
fc Levin] p005| , |Krumholz| p006| , KM06, and KMK08. 
Using simple models of core collapse in which angu- 
lar momentum is conserved in the collapse process and 
part of the matter is cast away by protostellar outflows 
( |Matzner fc McKee|2000| ), we find that both £ and T are 
higher in massive star formation than in low-mass star 
formation. In our models, the characteristic value of T 
rises from ~ 0.001 — 0.03 as one considers increasingly 
massive cores for which turbulence is a larger fraction of 
the initial support. 

The value of £ is more complicated, as it reflects the 
disk's thermal state as well as infalling acc retion rate, 
but the models of KMK08 and |Krumholz| ( |2006| indi- 
cate that its characteristic value increases from < 1 to 
~ 10 as one considers higher and higher mass cores - al- 
though the specific epoch in the core's accretion history 
is also important. In the case of ma ssive stars, such rapid 



accretion has been observed as in Beltran et al. ( 2006 1 



and Barnes et al. ( 2008 ) . Numerical simulations also find 



rapid accretion rates from cores to disks . Simulations 
such as those of Banerjee & Pudritz (20071 report £ ~ 10 
at early times in both magnetized and non-magnetized 
models. We note that T has significant fluctuations from 
core to core when turbulence is the source of rotation, 
and both £ and T are affected by variations o f the core 
accretion rate a round its characteristic value (Foster & 
Chevalier|p93 l. 

A major goal of this work is to probe the evolution of 
disks with £ > 1, as m ass accretion at this rate can not be 
accommodated by the Shakura & Sunyaev ( 1973 ) model 
with a < 1. Values of a exceeding unity imply very 
strong GI, and possibly fragmentation. 

3. NUMERICAL METHODOLOGY 

3.1. Numerical Code 

We use th e code ORION to conduct our numerical 
exper iments ( |Truelove et al. 1998 Klein 1999 Fisher 
|2002| . ORION is a parallel adaptive mesh refinement 
(AMR), multi-fluid, radiation- hydrodynamics code with 
self-gravity and lagrangian sink particles (Krumholz et 
al. 2004). Radiation transport and multi- fluids are not 
used in the present study. The gravito-hydrodynamic 
equations are solved using a conservative, Godunov 
scheme, which is second order accurate in both space 
and time. The gravito-hydrodynamic equations are: 

d : p= -V • (pv) - £ M t W{K - x,) 



d_ 
dt 



dt 



dt' 

i 

(pv) = -V • (pvv) - VP - pV(j) 
-^PiW(x-Xi) 

i 

(pe) = -V • [(pe + P)v] + pv • V(f> 



(3) 



(4) 



(5) 



^£ 4 iy(x-x 4 



Equations |3|-([6]) are the equations of mass, momen- 
tum and energy conservation respectively. In the equa- 
tions above, Mj, p^, and describe the rate at which 



mass and momentum are transfered from the gas onto 
the zth lagrangian sink particles. Summations in these 
equations are over all sink particles present in the cal- 
culation. W(x) is a weighting function that defines the 
spatial region over which the particles interact with gas. 
The corresponding evolution equations for sink particles 



dt 

~dt 
d_ 



M M, 

M 4 



-MjV0 - 



(G) 
(7) 
(8) 



These equations describe the motion of the point parti- 
cles under the influence of gravity while accreting mass 
and momentum from the surrounding gas. 

The Poisson equation is solved by multilevel elliptic 
solvers via the multigrid method. The potential 4> is 
given by the Poisson equation 



4ttG 



P 



^M^(x- Xi ) 



and the gas pressure P is given by 



P 



pfc^Tg 



(7-l)P 



e v 

2 



(9) 



(10) 



where T g is the gas temperature, p is the mean particle 
mass, and 7 is the ratio of specific heats in the gas. We 
adopt p — 2.33toh, which is appropriate for standard 
cosmic abundances of a gas of molecular hydrogen and 
helium. 

We use the sink par ticle implementation described in 
Krumholz et al. ( 2004 ) to replace cells which become too 
dense to resolve. Sink particle creation and AMR grid 
refinement are based on the Truelove criterion ( jTruelove] 
et al.|1997| which defines the maximum density that can 
be well resolved in a grid code as: 



P < Pj 



(11) 



G(Ax') 2 ' 

where, Nj is the Jeans number, here set to 0.125 for re- 
finement, and 0.25 for sink creation, and Ax' is the cell 
size on level /. When a cell violates the Jeans criterion, 
the local region is refined to the next highest grid level. 
If the violation occurs on the maximum level specified in 
the simulation, a sink particle is formed. Setting Nj to 
0.125 is also consistent with the resolution criterion in 



Nelson (2006). Sink particles within 4 cells of each other 
are merged in order to suppress unphysical n-body inter- 
actions due to limited resolution. At low resolution, un- 
physical sink particle formation and merging, can cause 
rapid advection of sink particles inwards onto the cen- 
tral star, generating spurious accretion. Moreover, be- 
cause an isotherm al, rotating gas filam ent will collapse 
infinitely to a line ( Truelove et al.|1997 l, an entire spiral 
arm can fragment and be merged into a single sink par- 
ticle. To alleviate this problem, we implement a small 
barotropic switch in the gas equation of state such that 

7=1.0001, p < p JS /4 (12) 
7=1-28, PJ s/A<p<pjs, (13) 
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where the J s subscript indicates the Jean's criterion used 
for sink formation. With this prescription, gas is almost 
exactly isothermal until fragmentation is imminent, at 
which point it stiffens somewhat. This modest stiffen- 
ing helps turn linear filaments into resolved spheres just 
prior to collapse and provides separation between new- 
born sink particles. The primary effect of this stiffening 
is to increase the resolution of the most unstable wave- 
length in a given simulation, at the expense of some dy- 
namical range. We describe the influence of this stiffen- 
ing on our results in £6.1 where we conduct some exper- 
iments in which it is turned off. 

As described via equations ([3])-([6]), sink particles both 
accrete from and interact with the gas and each other 
via gravity. Accretion rates are computed using a mod- 
ified Bondi-Hoyle formula which prevents gas which is 
not gra yitationally bound to th e p articles from accret- 
ing. See |Krumholz et al. | pOOl} and jOffner et al.| ( |2008 1 
for a detailed study of the effects of sink particle parame- 
ters. Note that we also use a secondary, spatial criterion 
for AMR refinement based on an analy tic p rediction for 
the disk size as a function of time (see [3.3|. 



3.2. Initial conditions 
We initialize each run with an isothermal core: 



p(r) 



Ac* 



(14) 



4ttGV 2 ' 

There is a small amount of rotational motion in our initial 
conditions, but no radial motion. A core with this profile 
is out of virial balance when A > 2, and accretes at a rate 

M_<corc f 0.975 (A = 2) 

M ~ G \ {2A)^/ir. (A » 2) (15) 

The value for A — 2 represents the |Shu ( |1977 1 inside- 
out collapse solution, whereas the limit A 3> 2 is derived 
assuming pressureless collapse of ea ch mass s hell. It is 



possible to predict M analytically (Shu 19771, but in 
practice we initialize our simulations with a range of val- 
ues A > 2 and measure M just outside the disk. Because 
our equation of state is isothermal up to densities well 
above the typical disk density (c Si( j — c StCorc ), MGj(? a core 
is equivalent to our parameter £. 

In order to set the value of our rotational parameter T 
and hold it fixed, we initialize our cores with a constant, 
subsonic rotational velocity: 

1/3 



ft = 



2Ac, 



(16) 



where w is the cylindrical radius. We arbitrarily choose a 
constant velocity rather than rigid rotation on spheres in 
order to concentrate accretion near the outer disk radii. 
Our definition of T in terms of the mean value of ji n 
rather than its maximum value is intended to reduce the 
sensitivity of our results to the choice of rotational pro- 
file. 

Given these initial conditions, our parameters £ and 
r remain constant throughout the simulation, while the 
collapsed mass and disk radius (as determined by the 
Keplerian circularization radius of the infalling material) 
increase linearly with time. We define a resolution pa- 
rameter, 

A = j^, (17) 



to quantify the influence of numerics on our results. Be- 
cause we hold the minimum grid spacing dx mm constant, 
A increases oc t as the simulation progresses. 

By artificially controlling the infall parameters of our 
disks, and then watching them evolve in resolution, we 
gain insight into the physical behavior of accretion with 
certain values of £ and T, as captured in a numerical 
simulation with a given dynamical range (A). Our initial 
conditions are necessarily ideal, allowing us to perform 
controlled experiments: of course realistic star-forming 
cores will undoubtedly be somewhat turbulent with time 
variable infall rates. 

3.3. Domain and Resolution 

Due to the dimensionless nature of these experiments, 
we do not use physical units to analyze our runs. The 
base computational grid is 128 3 cells, and for standard 
runs we use nine levels of refinement, with a factor of 
two increase in resolution per level: this gives an effec- 
tive resolution of 65,536 3 . More relevant to our results, 
however, is the resolution with which our disks are re- 
solved: A < 10 2 . To compare this to relevant scales in 
star formation, this is equivalent to sub-AU resolution in 
disks of - 50 - 100 AU. 

The initial core has a diameter equal to one half of the 
full grid on the base level. The gravity solver obeys pe- 
riodic boundary conditions on the largest scale; as the 
disk is 2.5 to 3 orders of magnitude smaller than the grid 
boundaries, disk dynamics are unaffected by this choice. 
The initial radius of the curr cnt infall is (vrr)- 2 / 3 i? kiin 
(from equations |2]), ( Jl4| , and (15l); although this is 
much larger than the disk itself, it is still ~ 15 — 40 
times smaller than the initial core and ~ 30 — 80 times 
smaller than the base grid. Tidal distortions of the in- 
fall are therefore very small, although they may be the 
dom inant seeds for the GI. We return to this issue in 
§6.3| where we compare two runs in which only the tidal 
effects should be different. 

In addition to the density criterion for grid refinement 
described in <j3] we also refine spatially to ensure that 
the entire disk is resolved at the highest grid level. We 
use £ and T to predict the outer disk radius (see $4} , and 
refine to 150% of this value in the plane of the disk!~~In the 
vertical direction, we refine to 40% of the disk radius: this 
value is larger than the expected scalehcight for any of 
our disks by at least ~ 15%. We find that we accurately 
capture the vertical and radial extent of the disk with 
this prescription, and the density criterion ensures that 
any matter at disk-densities extending beyond these radii 
will be automatically refined. 

3.4. Dynamical Self- Similarity 

Because our goal is to conduct a parameter study iso- 
lating the effects of our parameters £ and V, we hold 
each fixed during a single run. At a given resolution A, 
we expect the simulation to produce consistent results 
regarding the behavior of the accretion disk, the role of 
the GI, and fragmentation into binary or multiple stars. 
Since A increases linearly in time, each simulation serves 
as a resolution study in which numerical effects diminish 
in importance as the run progresses. Because the GI is an 
intrinsically unsteady phenomenon, a disk should fluctu- 
ate around its mean values even when all three of T, £, 
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and A are fixed. Because of this, and because A changes 
over the run, we expect our runs to be self-similar, but 
only in a limited, statistical sense. 

To illustrate how we expect resolution to affect our re- 
sults, consider another unsteady problem: the turbulent 
wake of a solid body in air. Our parameters £ and T are 
analogous to the macroscopic parameters of the problem, 
such as the body's aspect ratio and the Mach number of 
its motion, whereas A is analogous to the Reynolds num- 
ber of the flow. With all the parameters fixed, the flow 
field fluctuates around a well-defined average state. The 
average flow properties do depend on Reynolds number, 
but increasingly slowly in the high- Reynolds limit. In our 
simulations the resolution parameter is never fixed, but 
increases linearly in time; therefore we expect the disk to 
settle into an approximate steady state in which each ten 
orbits resemble the previous ten. Unlike a Sedov blast 
wave, we should not expect our disks to be exactly in- 
variant under scaling: this would not be consistent with 
the turbulent saturation of non-linear instabilities. 

Moreover, whereas many physical systems are captured 
perfectly in the limit of infinite resolution (A — > oo), this 
is not true of isothermal, gravitational gas dynamics, in 
which the mini mum mass and spacing o f fragments both 
scale as A -1 ( jlnutsuka fc Miyama|l992| . For this reason 
we quote the resolution A whenever reporting on the state 
of the disk-star system. 

We note that there exists a minimum scale in real ac- 
cretion disks as well, namely th e opacity-limited mini- 
mum fragment mass (| Rees|[l976| . The finite dynamical 
range of our numerical simulations is therefore analogous 
to a phenomenon of Nature, albeit for entirely different 
reasons. 

4. DISK PROPERTIES IN TERMS OF THE ACCRETION 
PARAMETERS 

To assess the physical importance of £ and T, it is use- 
ful to consider the case of a single star and its accretion 
disk. Because many £, T pairs lead to fragmentation, this 
assumption is only self-consistent within a subregion of 
our parameter space; nevertheless it helps to guide our 
interpretation of the numerical results. In order to as- 
sociate results from our parameters with those of previ- 
ous studies, we also derive expressions for disk averaged 
quantities such as Q and the disk-to-system mass ratio, 
fi as a function of £ and T. 

The combination 

\ _ 0)in C s,d _ Cs.d 

A) GAU d v Kin 1 ' 

is particularly useful, since it provides an estimate for 
the disk's aspect ratio (the scale height compared to the 
circularization radius). Being independent of M, it is 
more a property of the disk than of the accretion flow. 

The other important dimensionless quantity whose 
mean value depends primarily on £ and T (and slowly 
on resolution) is the disk-to-system mass ratio 

A* = T7 ■ (19) 

When the disk is the sole repository of angular mo- 
mentum, the specific angular momentum stored in the 
disk is related to the infalling angular momentum via: 

k = ( j^hr-) % (20) 



(j)in M *d 



where J; n is the total angular momentum accreted, so 
that Ji„/ ( (j) in M*d) = + 1) in an accretion scenario 



where 

(j) in /(2/i). 
can define 



cx MJ d . In our simulations U 



l, so jd 



Given the relation between jd and (j) in , we 



(21) 



Rd = [(lj + l)tA #k,in 
Vd = [(Z.+l^ftlcin 

which relate the disk's characteristic quantities (not the 
location of its outer edge) to conditions at the current cir- 
cularization radius i?k,in = (j)m /(GM»j). Such "charac- 
teristic" quantities are valuable for describing properties 
of the disk as a whole, rather than at single location, 
with an effective mass weighting. If we further sup- 
pose that the disk's column density varies with radius 
as S(r) cx r~ fcs (we expect fc E ~ 3/2 for a constant Q, 
isothermal disk) , we may define its characteristic column 
density E«j = (1 - k^/2)M d /(nR 2 d ): 



h 



(22) 



where/: 



(1 - fc E /2)(l + lj) 4 /n. Using equations (£8 1 
and (|21|)-(|22 1, we can rewrite the Toomre stability pa- 



rameter Q (ignoring the difference between ft and the 
epicyclic frequency for simplicity): 

Cc/v Cq^l d /~~\ 

" (23) 



Qa 
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(24) 



where Jq 
expect Qd 



(1 - fe E /2)(l + To the extent that we 
1 in any disk with a strong GI, this suggests 
fi ~ (r/£) 1/6 (l - fe E /2)- 1 / 2 (l + lj)- 1/2 ; and because we 
expect that \i has an upper limit of around 0.5 (see ^5] and 
discussion in KMK08 and |Shu et al.|1990 ), we see there is 
an upper limit to £ /T above which the system is likely to 
become binary or multiple. This is not surprising, as fj, is 
prop ortional to scale height when Q is constant; equation 
( 23 1 simply accounts self-consistently for the fact that fj, 
also affects R d . 

To go any f urther with analytica l argu ments, we must 
introduce the Shakura & Sunyaev ( 1973 1 a viscosity pa- 
rameterization, in which steady accretion occurs at a rate 

Mr) = ^f)r C -^ (25) 



Using the definition of £ 



Q(r) G 
3a(r) c s (r) 3 



CM < d (26) 

Insofar as Q ~ 1 when the GI is active, the effective value 
of a induced by a strong GI is directly proportional to £. 

The magnitude of T has important implications for disk 
evolution. As discussed previously by KMK08, V (called 
5i; n there) affects \x through the relation 



/' 



k.in 



Til 



1 



M, 



MM 



-1 -3(l-^)(l + I 3 > 



(27) 
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Fig. 1, — Two examples of single, binary, and multiple systems. The resolution across each panel is 328x328 grid cells. The single runs 
are £ = 2.9, T = 0.018 (top), £ = 1.6, T = 0.009 (bottom). The binaries are £ = 4.2, F = 0.014 (top), § = 23.4, T = 0.008, (bottom). The 
multiples are § = 3.0, T = 0.016 (top), £ = 2.4, V = 0.01 (bottom). Black circles with plus signs indicate the locations of sink particles. 
These correspond to runs 5, 1, 9, 16, 7, and 4 respectively. 



mass and angular momentum inflow, allowing the disk to 
drain and reducing its relative mass. We use the above 
relations to guide our interpretation of our simulation re- 
sults, specifically the dependence of disk parameters like 
/x, Qd, ct 7 and the fragmentation boundary, on £ and T. 

5. RESULTS 

Each of our simulations produces either a disk sur- 
rounding a single star, or binary or multiple star system 
formed via disk fragmentation; Figure [IJdepicts examples 
of each outcome. We use these three possible morpholo- 
gies to organize our description of the simulations. We 
explore the properties of each type of disk below as well 
as examine the conditions at the time of fragmentation. 

The division between single and fragmenting disks in 
£ and T is relatively clear from our simulation results, as 
shown in Figure [2j Several trends are easily identified. 
First, there is a critical £ beyond which disks fragment 
independent of the value of T. Below this critical £ value, 
there is a weak stabilizing effect of increasing T. As £ in- 
creases, disks transition from singles in to multiples, and 
finally into binaries. We discuss the distinction between 
binaries and multiples in |5.3| 

In table [l] we list properties of the final state for all of 
our runs, their final multiplicity (S, B, or M for single, bi- 
nary, or multiple, respectively), and the disk-to-star mass 
ratio /i / measured at the time at which we stop each ex- 
periment, as well as the maximum resolution A„. Note 



where the second line uses disk-averaged quan titie s to 
construct a mean accretion rate from equation ( 25 ) . In 
our simulations /i ~ so we expect to saturate at the 
valu e for which the two terms on the right of equation 
(27) are equal, 

^V? + 2B)V»-B, where B = ^j^g^ - 

(28) 

The disk mass fraction increases with B, so both T 
and £ have a positive effect on /i, whereas a tends to 
suppress the disk mass. Note that, when B is small and 
/i ~ V%B, equation (23 1 implies Qd ~ 3a/£ in accor- 



dance with equation (25). Because the effective value of 



a induced by the GI is a function of disk parameters, we 
cannot say more without invoking a model for a(T, £) or 
a(Q,n) as in KMK08. 

The scalings of disk properties with the dimensionless 
parameters of the problem are in accord with intuitive ex- 
pectations. An increase in £ corresponds to an increase 
in accretion rate at fixed disk sound speed, and as a re- 
sult the equilibrium disk mass rises. An increase in T 
corresponds to an increase in the mean angular momen- 
tum of the infall at fixed sound speed, leading to larger 
disks that must transport more angular momentum, and 
thus again become more massive. An increase in a cor- 
responds to an increase in the rate at which the disk can 
transport angular momentum and mass at a fixed rate of 
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that the disk extends somewhat beyond i?k,in: therefore 
the disk as a whole is somewhat better resolved than the 
value of X n would suggest. For the disks which fragment, 
we also list the value of /x/, A/ and Q just before frag- 
mentation occurs. 
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TABLE 1 



Each run is labelled by £,r, multiplicity outcome, the 

FINAL VALUE OF THE DISK-TO- STAR MASS RATIO ,/i AND THE FINAL 
RESOLUTION, A n . VALUES OF V ARE QUOTED IN UNITS OF 10~ 2 . 
FOR FRAGMENTING RUNS THE DISK RESOLUTION Xf, Q 2 D 
(EQUATION |20p AND jij AT THE TIME OF FRAGMENTATION ARE 

LISTED AS WELL. S RUNS ARE SINGLE OBJECTS WITH NO PHYSICAL 
FRAGMENTATION. B'S ARE BINARIES WHICH FORM TWO DISTINCT 
OBJECTS EACH WITH A DISK, AND M ARE THOSE WITH THREE OR 

MORE STARS WHICH SURVIVE FOR MANY ORBITS. * INDICATES RUNS 
WHICH ARE NOT SUFFICIENTLY WELL RESOLVED AT THE TIME OF 

FRAGMENTATION TO MAKE MEANINGFUL MEASURES OF fl f , AND Q. 



In table [2] we describe those disks which do not frag- 
ment: we list the analytic estimate for the characteristic 
value Too mrc 's Q, Qd, the measured minimum of Qid 
(equation |29| , the radial power law fcs which character- 
izes S(r) for a range of radii extending from the accre- 
tion zone of the inner sink particle to the circularization 



radius i?k,in, 
acteristic disk radius 



the final disk resolutio n, \ n , and the char- 
Rd (equation (21 1. 
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TABLE 2 

Single runs (numbers as from table [!}. We list values for 

the characteristic predicted value of toomre's q, q d 
(equation |23ll, a s well as the measured disk minimum, q2d 

EQUATI0n"1[29[| . WE ALSO LIST THE SLOPE OF THE SURFACE 
DENSITY PROFILE, fe s AVERAGED OVER SEVERAL DISK ORBITS, THE 
FINAL RESOLUTIONS, AND R d AT THE END OF THE RUN (EQUATION 

I2B 



5.1. The Fragmentation Boundary and Q 

It is difficult to measure a single value of Q to char- 
acterize a disk strongly perturbed by GI, so we consider 
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1 10 

Fig. 2. — Distribution of runs in £ — T parameter space. The 
single stars are confined to the low £ region of parameters space, 
although increasing T has a small stabilizing effect near the transi- 
tion around £ = 2 due to the increasing ability of the disk to store 
mass at higher values of T. The dotted line shows the division be- 
tween single and fragmenting disks: V = £ 2,5 /850. As £ increases 
disks fragment to form multiple systems. At even higher values 
of £ disks fragment to make binaries. We discuss the distinction 
between different types of multiples in §5.3| The shaded region of 
parameter space shows where isothermal cores no longer collapse 
due to the extra support from rotation. 

two estimates: a two dimensional measurement Q2D, and 
a one-dimensional measure Q av (r) based on azimuthally- 
averaged quantities. 



?2D(r, 4>) : 



ttGE 1 
c s (r)R(r) 



(29) 
(30) 



(bars represent azimuthal averages). As Figure [3] shows, 
the two-dimensional estimate shows a great deal of struc- 
ture which is not captured by the azimuthal average, 
let alone by Qd- Moreover, while the minimum of the 
averaged quantity is close to two, the two dimensional 
quantity drops to Q ~ 0.3. We find that the best pre- 
dictor of fragmentation is the minimum of a smoothed 
version of the two-dimensional quantity (smoothed over 
a local Jeans length to exclude meaningless fluctuations) , 
although Qd shows a similar trend. We use this quantity 
in table [T] and compare it to the analytic estimate Qd in 
table [2] for non-fragmenting disks. 

We emphasize that the critical values of Q at which 
fragmentation sets in depend on the exact method used 
for calculation (e.g. Q az or Q2d)- Moreover, we do not 
expect to reproduce fragmentation at the canonical order 
unity boundary. This only marks the crit ical case for th e 
m=0 unstable mode in razor-thin disks ( |Toomre | |1964[ ). 
As discussed by numerous authors, the frag mentation 
criterion is somewhat different for thick disks ([G ol dreich| 
fc Lynden Bell|1965||Laughlin et al. |1997[ [1998D , and the 
growth of higher order azimutha l modes (Adams et al. 
1989||Shu et al.||1990| [Laughlin fc Korchagin||1996| ) 



Another consequence of trying to describe thick disks 
with multiple unstable modes is that the fragmentation 
boundary cannot be drawn in Q-space alone. We use 
Q2D and /1 in Figure [4] to demarcate the fragmentation 
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Fig. 3.— Top: Q av in a disk with £ = 2.9, T = 0.018. The 
current disk radius, -Rk.m is shown as well. Bottom: Log(Q 2 l?) 
(equation |29| l in the same disk. While the azimuthally averaged 
quantity changes only moderately over the extent of the disk, the 
full two-dimensional quantity varies widely at a given radius. Q is 
calculated using n derived from the gravitational potential, which 
generates the artifacts observed at the edges of the disk. Here and 
in all figures, we use 8 X to signify the resolution. 



boundary. Labeled curves illustrate that the critical Q 
for fragmentation depends on the disk scale height (equa- 



tion 18 1. At a given value of Q, a disk with a larger value 
of fi will have a larger aspect ratio, a nd w ill therefore be 
more stable. Recall from equation (18 1, that the disk 



aspect ratio is proportional to (^/r) 1 . 

This trend i s cons istent with the results of lGoldreich fcl 
Lynden-Bell (1965) for thick disks; because the column 
of material is spread out over a larger distance, H, its 
self-gravity is somewhat diluted. The fact that two pa- 
rameters are necessary to describe fragmentation is also 
apparent in Figure [2j where the boundary between sin- 
gle and multiple systems is a diagonal line through the 
parameter space. 

Although two criteria are necessary to prescribe the 
fragmentation boundary, we observe a direct correspon- 
dence between \i and T, and £ and Toomre's Q. Figure [5] 
shows that fi m 2T 1 / 3 for both single star disks, and just 
prior to the onset of fragmentation in disks that form 
binaries and multiples. We find a similar correspondence 
between £ and the combination Qd^, which is a direct 
correlation between £ and Q defined with respect to the 
disk circularization radius (using Rd in the definition of 
Qd brings in an extra factor of //.) 
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Q 

Fig. 4. — Steady-state and pre-fragmentation values of Q and 
/i for single stars and fragmenting dis ks respectively. We use the 
minimum of Q2D as described in §5,1| Symbols indicate the mor- 
phological outcome. Note that the non-fragmenting disks (large 
triangles) have the highest value of p, for a given Q. Contours 
show the predicted scaleheight as a function of Q and /i. It is clear 
that the single disks lie at systematically higher scale heights. We 
have assumed = 3/2 in calculating scaleheight contours as a 
function of Q and fj,. 

5.2. Properties of non-fragmenting disks 

Although we quote a single power law value for the sur- 
face density profiles of disks in table|2j the surface density 
structure is somewhat more complex. We find that the 
disks show some evidence of a broken power law struc- 
ture: an inner region, characterized by fcs, where disk 
material is being accreted inwards, and an outer region 
characterized by a steep, variable power law due to the 
outward spread of low-density, high angular momentum 
material. We find disks characterized by slopes between 
fcs = 1 — 2. Clustering around fcs = 3/2 is expected, as 
this is the steady-state slope for a constant Q, isothermal 
disk. Our measurements of Q(r) (equation 29 1 show fluc- 
tuating, but roughly constant value over the disk radius. 
Note that the slope of the inner disk region tends to in- 
crease with r. Figure [6] shows normalized radial profiles 
for the non- fragmenting disks. Profiles are averaged over 
approximately three disk orbital periods. The flattening 
at small radii is due to the increasing numerical viscosity 
in this region (j |6.3[ ). 

We find an upper mass limit of /1 ~ 0.55, for single 
stars, which means that disks do not grow more massive 
than their ce ntral star. A max imum disk mass has been 
predicted by |Shu et ah] ( 1990 1 as a consequence of the 
SLING mechanism. Such an upper limit is expected as 
eccentric gravitational instabilities in massive disks shift 
the center of mass of the system away from the central 
object. Indeed, we observe this wobble in binary forming 
runs. The subsequent orbital motion of the primary ob- 
ject acts as an indirect potential exciting strong m = 1 
mode perturbati ons which can induce binary formation 
(Shu et al. 1990). We find that this maximum value is 
consistent with their prediction. 

Using the analytic expressions above, we can also de- 
rive an expression for an effective Shakura-Sunyaev a. In 
this regime of parameter space, £ and T are always such 
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Fig. 5. — At right r ys (i with the fit in equation (|3ll overplotted. At left, Qd^ vs £ with the scaling Q <x £ 1 / 3 overplotted. Runs, 16, 
17, 18 are omitted as the low resolution at the time ot fragmentation makes measurements of fi and therefore unreliable. 



that B <C 1 (assuming a does not stray far from unity). 
We therefore expect that fi oc r 1 / 6 ^ 1 / 3 a~ 1 / 2 . Using this 
relation we can find a functional form of a(£,r). Our fit 
to the data shown in Figure [5] implies 



2r l/3 



(31) 



with some scatter for both single disks and fragmenting 
disks just prior to fragmentation. We can use th is fit to 
infer a scaling relation for a using equation ( 28 1 in the 
limit [i ~ y/2B: 



1 



ad 



2/3 



18(2-fc E ) 2 (l + Z,) 2 r 1 /3' 



(32) 



The scaling is consistent with our expectation that 
driving the disk with a higher £ causes it to process ma- 
terially more rapidly, while increasing T decreases the 
efficiency with which the disk accretes. Equation (32 1 
predicts disk averaged values of a for single star disks 
between ~ 0.3 — 0.8. These values are consistent with 
the obser ved accretion rates, and numerically calculated 
torques (5 5.4 1. 



5.3. The formation of binaries and multiples 

As shown by Figure 2] a large swath of our parameter 
space is characterized by binary and multiple formation. 
We find that the division between fragmenting and non- 
fragmenting disks can be characterized by a minimum 
value of r at which disks of a given £ are stable. In 
Figure[2]we have plotted this boundary as T — £ 25 /850. 

While we do not claim that our numerical experiments 
are a true representation of the binary formation process, 
we do expect to find binaries in much of the parameter 
space characteristic of star formation, as nearly half of all 



stars are in binaries (Duqucnnoy & Mayor 1991 Duchenc 



et al.|2007 1. Moreover, as the binary forming parameters 



are typical of higher mass star formation, where bina- 
ries and multiples are expected to comprise perhaps 75% 



of systems, these findings are encouraging ( |Mason et al. 
|1998| |2009[ ) . We discuss several general trends here, but 
defer a detailed analysis of binary evolution and applica- 
tion to observations to a later paper. 
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Fig. 6. — Normalized density profiles for the single-star disks. 
Profiles are azimuthal averages of surface densities over the final 
~ 3 disk orbital periods. We find that while the inner regions are 
reasonably approximated by power law slopes, the slope steepens 
towards the disk edge. For comparison, slopes of fcg = 1, 1.5, and 
2 are plotted as well. Runs are labelled according to their values 
in table HI 

Are these equal mass binaries? Low mass stellar com- 
panions? Or maybe even massive planets? In a self- 
similar picture it is difficult to tell. In an actively accret- 
ing multiple system, as long as the mass reservoir has an- 
gular momentum such that the circularization radius of 
the infalling material is comparable to the separation be- 
tween objects, the smaller object, which is further from 
the cent er of mass, will accrete due to the torque im - 



balance (Bate & Bonnell 1997 
Similarly, in thick, gravitational 



lation mass approaches the stellar mass: 



Bonnell fc Bate||1994b I 
y unstable disks, the iso- 



M iso = 4Trf H r H r d j: « 30 



3.5 



3/2 



Q~ 3/2 M*. (33) 
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Here, r H = (M S /3AQ 1/3 is the Hill radius, M s , and M* 
are the masses of the secondary and primary, and the 
numerical factor fu, represents how many Hill radii an 
object can feed fro m in the disk - numerical si mulations 
suggest f H ~ 3.5 flLissauer||T987l |Rafikov||2002| ). There- 
fore the evolution of these objects in our models is clear: 
they tend to equalize in mass. The binary separation 
will also grow if any of the infalling angular momentum 
is transferred to the orbits as opposed to the circumstel- 
lar disks. These trends are borne out in our experiments: 
binary mass ratios asymptote to values of 0.8 — 0.9 and 
separations to ~ 60% of i?k,in- In a realistic model for 
star formation, the parameters that characterize a single 
run in this paper will represent only one phase in the life 
of a newborn system. The trajectory through £ — T space 
which the systems take following binary formation will 
strongly influence the outcome in terms of separation and 
mass ratio. For example, should the disk stabilize and 
accretion trail off quickly following binary formation, it 
is quite likely that a large mass ratio would persist as the 
disk drains preferentially onto the primary object once 
the secondary reaches its isolation mass. By contrast, in 
systems which fragment before most of the final system 
mass has accreted, we expect more equal mass ratios. 

5.3.1. Hierarchical multiples and resolution dependence 

Disks which are at the low £ end of the binary forming 
regime tend to form binaries at later times, and therefore 
at higher disk resolution. One consequence of this is the 
formation of hierarchical multiples. When disks become 
violently unstable, they fragment into multiple objects. 
Because of the numerical algorithm which forces sink 
particles within a gravitational softening length of each 
other to merge, at lower resolution many of these parti- 
cles merge, leaving only two distinct objects behind. At 
higher resolution, while some of the particles ultimately 
merge, we find that three or four objects typically survive 
this process. We cannot distinguish between merging and 
the formation of very tight binaries. In addition to merg- 
ing, small mass fragments are occasionally ejected from 
the system entirely. This appears be a stochastic pro- 
cess, though we have not done sufficient runs to confirm 
this conclusion. 

Disks which form binaries at early times and develop 
two distinct disks can also evolve into multiples when 
each disk becomes large enough and sufficiently unstable 
to fragmentation. In general, once a binary forms, the 
system becomes characterized by new values of £ and T 
which are less than those in the original disk. As the 
distribution of mass and angular momentum evolves in 
the new system, the relative values of £ and T evolve as 
well. However, once the mass ratios have reached equi- 
librium as is the case for run #16 shown in the bottom 
center of Figure [TJ each disk sees £ of roughly half the 
original value, which for an initial £ ~ 24 is still well 
into the fragmenting regime. As a result, the fact that 
the two disks ultimately fragment is expected. On the 
contrary, for the lowest £ binary runs, once one fragmen- 
tation event occurs, the new £ may be sufficiently low 
to suppress further fragmentation. The evolution of T in 
the newly formed disks is more complicated, depending 
on how much angular momentum is absorbed into the 
orbit as compared to the circumstellar disks. We defer a 
discussion of this to a later paper. 
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Fig. 7. — Azimuthal averages of d iffer ent components of torque 
expressed as an eff ecti ve a (equation |35| for run #8. The straight 
line, ay (equation \'>2\ is plotted for comparison. The agreement 
between the analytic value of ay and the combined contribution 
from the other components is best near the expected disk radius 

It is clear that there is a numerical dependence to this 
phenomenon which we discuss in |6.3[ but there is a cor- 
respondence with the physical behavior of disks as well. 
The radius and mass of a fragmenting disk are likely 
to influence the multiplicity outcome of a real system. 
Cores with high values of £ that form binaries early in our 
numerical experiments correspond to cores whose disks 
fragment into binaries at small physical size scales, where 
the disk may only be a few fragment Hill radii wide, and 
contain a relatively small number of Jeans masses. It is 
possible that at these size scales, numerous bound clumps 
in a disk might well merge leaving behind a lower multi- 
plicity system than one in which the ratio of disk size to 
Hill radii or mass to Jeans mass is higher. 

5.4. Gravitational Torques and Effective a 

We verify that the accretion observed in our disks 
is generated by physical torques by computing the net 
torque in the disk. It is convenient to analyze the torques 
in terms of the stress tensor, Tr^, which is made up of 
two components: large scale gravitational torques and 
Reynolds stresses. Following Lodato & Rice (2005 1 we 
define: 



Tr4 — 



9R9± 
4ttG 



(34) 



where 5v — v — v. In practice, we set Jvr = ~vr, while 
(5v0 is calculated with respect to the azimuthal average of 
the rotational velocity at each radius. In reality there is 
an extra viscous term attributable to numerical diffusion. 



We discuss the importance of this term in {6.3 



The first term in equation ( 34 1 represents torques due 
to large scale density fluctuations in spiral arms, while 
the second is due to Reynolds stresses from deviations 
in the velocity field from a Keplerian (or at least radial) 
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velocity profile. To facilitate comparison with analytic 
models, the torques can be represented as an effective a 
where: 

dlnfi 



T, 



R4> 



dlni? 



(35) 



We can compare these torques to the characteristic disk 
a c i in equation ( 32 1 . Although there is variability in the 
disk accretion with time, it is consistent with a constant 
rate over long timescales. 

Figure [7] compares ay to the azimuthal average of the 
physical torques for one of our runs. We also show t he ex - 
pected contribution from numerical diffusion (see £6.31. 
The accretion expected from these three components is 
consistent with the time averaged total accretion rate 
onto the star. Due to the short term variability of the 
accretion rate, the two do not match up exactly. It is in- 
teresting to note the radial dependence of the Reynolds 
stress term, which in the inner region decays rapidly, 
before rising again, due to the presence of spiral arms. 
In both the azimuthal average and the two dimensional 
distribution we see that at small radii numerical diffu- 
sion dominates, whereas at large radii deviations in the 
azimuthal velocity which generate Reynolds stresses are 
spatially correlated with the spiral arms. 

5.5. Vertical Structure 

When the disks reach sufficient resolution, we can re- 
solve the vertical motions and structure of the disk. We 
defer a detailed analysis of the vertical structures to a 
later paper, but discuss several general trends here. De- 
pending on the run parameters, the disk scale height is 
ultimately resolved by 10-25 grid cells. We observe only 
moderate transonic motions in the vertical direction of 
order M. ~ 1 — 2. Figure [8] shows two slices of the z- 
componcnt of the velocity field for a single system, one 
through the X-Z plane, and the other through the disk 
midplanc. Although there is significant substructure, the 
motions are mostly transonic. 

We also observe a dichotomy in the vertical structure 
between single and binary disks. Although the values 
of £ and T should dictate the scaleheight (see equation 
(19}), and therefore higher £ disks which become binaries 
should have smaller scalcheights to begin with, we ob- 
serve a transition in scaleheight when a disk fragments 
and becomes a binary. Large plumes seen in single disks, 
like those shown in Figure [9] contain relatively low den- 
sity, high angular momentum material being flung off 
of the disk. The relatively sharp outer edges are cre- 
ated by the accretion shock of infalling material onto 
these plumes. We observe small scale circulation patterns 
which support these long lived structures. Disks sur- 
rounding binaries, by comparison remain relatively thin; 
in particular while the circumprimary disks are slightly 
puffier than expected from pure thermal support, the cir- 
cumbinary disk (when present) is sufficiently thin that we 
do not consider it well resolved. This implies that the ef- 
fective r values that binary disks see declines more than 
£ according to equation (18). This is consistent with the 
statement that some of the infalling angular momentum 
is transferred into the orbit instead of on to the disks 
themselves. 

6. CAVEATS AND NUMERICAL EFFECTS 
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Fig. 8. — Cuts along the vertical axis and disk midplane of the 
vertical velocity, normalized to the disk sound speed. Clearly most 
of the vertical motions in the disk are transonic, although at the 
edges of the disk the velocities exceed M ~ 1. 

6.1. Isothermal equation of state 

Many simulations have shown the dra matic effects that 
thermodynam i cs have on disk behavior fBoss et al.|2000| 
Gammie|200TI |Rice et al.|20051 ILodato fe Rice|2005| |ljo-| 
ley et al .|20"06] |Krumholz et al.|2007||Offner et al.|2009| . 

Since we are concerned with fragmentation, we must be 
aware of the potential dependencies of the fr agmenta- 
tion b oundary on cooling physics. Starting with |Gammie| 
(20011, there has been much discussion of the "cooling 
time constraint" that states that a disk with Q ~ 1 will 
only fragment if the cooling time is short. While this is 
a valuable analysis tool for predicting the evolution of 
a system from a snapshot and for quantifying the feed- 
back from gravito-turbulence, for most of the protostellar 
disks that we are modeling, the cooling time at the lo- 
cation of fragmentation is short. In the outer radii of 
protostellar disks in general, irradiation is the dominant 



source of heating dD'Alessio et al.|1997||Matzner fc Levin 
2005| |Krumholz et al.|2007[ |Kratter et al.|2008p . In fact, 
even a low temperature radiation bath c an contribute sig - 



nificantly to the heat budget of disks ( Cai et al7|[2008 1 
Such passively heated (through irradiation) disks behave 
more like isothermal disks than barotropic disks, because 
the energy generation due to viscous dissipation is small 
compared to the energy density due to radiation. Con- 
sequently, feedback from accretion in the midplane does 
not alter the disk temperature significantly. For realistic 
opacity laws, disks which are dominated by irradiation 
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cool on timescales much faster than the orbital period 
at large radi i. Numerical simulations such as |Krumholz| 
et al. (2007) find that strongly irradiated disks have a 
nearly isothermal equation of state. In fact, the mor- 
pholog ical outcome is similar to those of |Krumholz et ah] 
( |2007| with comparable values of £. 

Another possible concern is the lack of a radial tem- 
perature gradient, independent of the equation of state. 
Both passively and actively (through viscous dissipation) 
heated disks will be warmer at small radii, though the 
radial dependency changes with the heating mechanism. 
Actively heated disks typically have steeper gradients. 
In cither case, it is possible that the warmer inner disk 
would be stabilized against gravitational instability and 
slow down accretion. In these experiments, we find spiral 
arms persist in regions where the average value of Q is 
well above that at which instability is presumed to set 
in, with local values exceeding this by an order of mag- 
nitude. It seem plausible that due to the global nature 
of the low- to spiral modes, angular momentum trans- 
port may still occur in regi ons one would assum e stable 
against GI. As discussed by Adams et al. (19891, to = 1 
modes can have appreciable growth rates for remarkably 
high values of Q when the evanscent region is as much 
as 70% of the disk radius. In the event that the GI does 
shut off due to increasing temperature, material from 
the outer, unstable portion of the disk will likely accu- 
mulate until the critical surface density for GI is reached, 
causing ky, to steepen. Further numerical investigation 
of this is necessary; we point readers to high resolution 
st udies of disks with rad i al temperatu r e gradients such 



Krumholz et al. (2009); Boleyet al. (20071, Cai et al 



as ivrumnolz et al. (zuuy); fiolcy et 
( [2008P , and |Krumholz etaL| |2007p . 
In order to test the effects of the g£ 



gas stiffening we intro- 
duced to avoid unphysical merging of our sink particles 
(see (3.1|, we have conducted several purely isothermal 
experiments in which it is turned off. The removal of 
the barotropic switch aritifically enhances accretion at 
early times due to sink particles formed via numerical 
fragmentation merging with the central star. Remov- 
ing the barotropic switch is equivalent to increasing the 
resolution of the fragmentation process, but decreasing 
the resolution of the scale of fragmentation relative to A, 
the disk resolution. Using a barotropic switch allows the 
disk to reach a higher A before fragmentation sets in for 
a given set of parameters. 

6.2. Insensitivity of disk dynamics to core temperature 

Our parameterization of disk dynamics is based on the 
idea that thermodynamics can be accounted for by one 
parameter, £, which compares the accretion rate to the 
disk sound speed c s ^. A basic corollary of this notion 
is that the core temperature c s . cor c has no effect on disk 
dynamics, except insofar as it affects the accretion rate. 
We have defined £ with respect to the disk sound speed, 
but since our disks and cores are the same temperature, 
we could equally well have used the core sound speed. 
Therefore the question arises whether £ should be com- 
puted by normalizing the accretion rate to the disk sound 
speed, c s ,d, or the core sound speed, c S;C . To test this, 
we ran simulations in which c St d and c s . core differed: we 
imposed a change in temperature over a range of radii in 
which infall is highly supersonic. 

To demonstrate that £ defined with respect to the disk 
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Fig. 9. — Density slices showing vertical structure in a single and 
binary disk. The top plot is a single star with £ = 1.6, T = 0.09, 
while the bottom is a fragmenting binary system with £ = 24.3, T = 
0.008. The extended material in the binary system is generated by 
a combination of large scale circumbinary torques and the infalling 
material. Colorscale is logarithmic. The box sizes are scaled to 
1.5-Rk,in hi the plane of the disk. 

sound speed is indeed a better predictor of the morpho- 
logical and physical behavior of the disk, we compare 
A at the time of fragmentation, A/ to both £ and the 
equivalent parameter defined in terms of the core sound 



speed, £ c 



GM/ci 



for runs with similar values of 



r. We observe a correlation between resolution at the 
time of fragmentation and £ at fixed T, and so if core 
temperature is irrelevant, these runs should follow the 
same trend. 

Figure [To] shows that A/ correlates extremely well with 
£ at similar values of T, but poorly with £ cor c for the 
heated runs. The scaling of Xf with £ is also related 
to the existence of an upper limit on fj, as a function of 
T: disks with higher £ approach this critical value of /x 
faster, and thus at lower A. 

6.3. Resolution 



We have shown in §5.4| that the observed accretion is 
consistent with the combined gravitational torques and 
Reynold stresses, and that these are dominant over that 
expected purely from numerical diffusion. Because of the 
self-similar infall, convergence to a steady state within 
a given run is a good indicator that numerics are not 
determining our result; in effect, every run is a resolu- 
tion study. That we observe a range of behavior at the 
same resolution but different input parameters also im- 
plies that numerical effects are sub-dominant. We con- 
sider our disks to begin to be resolved when they reach 
radii such that R^^/Ax > 30. The effective numerical 
dif fusivity, which we plot i n Figure [7] has been estimated 
by |Krumnolz et al.|p004| for ORION. Specifically they 
find that: 

tb (^r 85 . (36) 



78 



Ax VAa 



where 



GAL 



TB = 



(37) 



is the standard Bondi radius. 

For our typical star and disk parameters, this implies 
numerical a's of order 0.1 — 0.3 at the minimum radius 
at which we are resolved. This implies that for our "low" 
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Fig. 10. — Correlation between \j and the infalling accretion rate 
for heated and non heated runs with comparable V. Plus symbols 
indicate non-heated runs, and the crosses are heated runs. The 
arrows and red crosses indicate the position of the runs evaluated 
with respect to § coro . Runs shown have V values ranging from 0.006 
to 0.009. The shaded region illustrates the scaling \j oc £ — x . This 
scaling is related not only to the existence of a critical value of fi, 
but also tied to the effect of resolution on fragmentation. 

accretion rate cases, at most 1/3 of our effective alpha 
could be attributed to numerical effects at l ow resolu- 
tion. See discussions by Offner et al. (2008); Krumholz 



et al. (2007 2004) for a detailed analysis of disk reso 



iution requirements. At our resolution of 50-100 radial 
cells across the disk, the dominant effect of numerical 



diffusion is likely a suppression of fragmentation phetty 

177 



& Ostriker 2006| Nelson 20061. Because the isothermal 
spiral arms can become very narrow prior to fragmenta- 
tion, numerical diffusion across an arm may smear out 
some overdensities faster than they collapse. Therefore 
the conclusions regarding the fragmentation boundary 
are likely conservative. 

We explicitly demonstrate morphology convergence in 
one of our binary runs. We rerun run #16 (as labelled in 
table [TJ at double the resolution (128 3 with 10 levels of 
refinement as opposed to 9) . Increasing the physical res- 
olution also decreases the code time step proportionally 
so that the ratio of the timestep to orbital period as a 
function of A should be preserved. In fact, there is little 
that can be different between the runs at two resolutions 
at the same effective A. 

The two runs have the same morphology, and charac- 
teristic disks properties as a function of A, as expected. 
We show in Figure [IT] snaphots of the standard and high 
resolution runs. The standard resolution run (left) is 
at twice the elapsed "time" of the high resolution one 
(right), and so the same numerical resolution, A. We con- 
firm that the mass accretion rate is consistent between 
the two runs: at the snapshots shown the mass ratio of 
the lower resolution run is 0.46, while the higher resolu- 
tion run is 0.48. We consider this variation to be w ithin 
the expected variation of the parameters (see ^3.4 1. To 
the extent that numerical artifacts are seeding instabil- 



ities, we expect some stochasticity in the details of the 
fragmentation between any two runs. Although the ef- 
fect is small, it is also possible that since the physical size 
of the disk (and the radius from which material is cur- 
rently accreting) relative to the box size is larger at the 
same value of A for the low resolution run, the large scale 
quadrapole potential from the image masses is stronger 
in the low resolution case. 

We also compare a multiple run at a lower resolution 
by a factor of two. Again we find the same morphological 
outcome. We find that the disks behave equivalently at 
the same radial resolution although the elapsed time, and 
dimensional masses are different. 



The scaling of Ay with £ in Figure 10 also demonstrates 
that resolution plays a role in determining when disks 
fragment. Although the infall is self-similar, the disks ap- 
proach a steady state as parameters like Q and /i evolve 
toward constant values. This evolution, and sometimes 
fragmentation, is influenced by the interplay between de- 
caying numerical viscosity and increasing gravitational 
instability in the disk as a function of A. 

7. COMPARISON TO PREVIOUS STUDIES 

The literature is replete with useful simulations of 
protostellar and protoplanetary disks at various stages 
of evolution, however most involve isolated disks, with- 



out infall at large radii (|Laughlin fc Bodenheimer] 1994 
Laughlin & Rozyczka 1996; Rice et al. 2005, Lodato & 



Ricel|2 005 Fr omang fc Nelson||2006| |Shetty fe Ostriker 
2006] |Boley et al.||2006| |Lodato et al.||2007| |Cai et al. 

2008 ) . These simulations include a wide range of physics, 
from magnetic fields to radiative transfer, but due to the 
lack of infalling matter, they neither develop disk pro- 
files (surface density, temperature) self-consistenly, nor 
do they enter the regime of interest in this work: rapid 
accretion in the embedded phase. For a review of many 
of the issues addressed by current GI disk simulations, 
see|Durisen et al. (20071. 

here are a few simulations of self-consist ent growth 



T 



and evolution (iVorobyov & Basu||2007 20081. These are 



ideal for following the long term evolution of more qui 
escient lower mass disks. However, because they are two 
dimensional, and lack a moving central potential, they 
cannot follow the evolution of non-axisymmetric modes 
which are driven by the displacement of the central star 
from the center of mass, nor can they accurately simulate 
the formation of multiple systems. Other authors have 
investigated the initial stage of core collapse onto disk s 
( TBanerjee fc Pudritz|[2007] |Tsuribe fc Inutsuka|[l999a] , 
however these authors focus on the effects of magnetic 
fields and fragmentat ion of the core prior to disk for- 
mation respectively. Tsuribe fc Inutsuka ( |1999b| ) and 



Matsumoto fc Hanawa ( |2003 ) have also investigated the 
collapse of cores into disks and binaries, thou gh they do 

for detailed 

2007T 



not investigate many disk propert i es (see { 7. 1 
comparisons ) 



et al. (2009 



< 2 
fci 



Krumholz et al 
have conducted three 



Krumholz 



imcnsional radiative 



transfer calculations, but due to computational cost can 
only investigate a small number of initial conditions. 

In addition to numerical work, there are a range of 
semi-analytic model s which follow the time evolution 
of accreting disks (IHueso fc Guillot| [20051 KMK08). 
KMK08 examined the evolution of embedded, massive 
disks in order to predict regimes in which gravitational 
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Fig. 11. — At left: a snapshot of the standard resolution of run #16 shortly after binary formation. At right, the same run at double the 
resolution. Because of the self-similar infall prescription, we show the runs at the same numerical resolution, as time and resolution are 
interchangeable. In this case the high-resolution run has taken twice the elapsed "time" to reach this state. The two runs are morphologically 
similar and share expected disk properties. 

instability, fragmentation of the disk, and binary forma- 
tion were likely. They concluded that disks around stars 
greater than 1 — 2Af Q were likely subject to strong grav- 
itational instability, and that a large fraction of O and 
B stars might be in disk-born binary systems. |Hueso fc] 
Guillot (20051 have also made detailed models of disk 
evolution, though they examine less massive disks, and 
do not include explicitly gravitational instability, and 
disk irradiation. 

In KMK08 we hypothesized that the disk fragmenta- 
tion boundary could be drawn in Q- \i parameter space, 
where small scale fragmentation was characterized by low 
values of Q and binary formation by high values of fi. 
Due to the self-similar nature of these simulations, the 
distinction between these two types of fragmentation is 
difficult, as the continued accretion of high angular mo- 
mentum material causes the newly formed frag ment to 
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preferentially a ccrete material and grow in mass (Bonnell 
fc Bate|1994a [. Moreover, because the disks are massive 



and thick, the isolation mass of fragments is compara- 
ble to the disk mass, and so there is little to limit the 
continued growth of fragments. 

7.1. The evolution of the accretion parameters in the 
isothermal collapse of a Bonnor-Ebert Sphere 

While self-similar scenarios are useful for numerical ex- 
periments, they do not accurately capture the complex- 
ities of star formation. In particular, in realistic cores, 
£ and T evolve in time. Therefore it is interesting to 
chart the evolution of a more realistic (though still ideal- 
ized) core through our parameter space. We consider the 
isothermal collapse of a Bonnor-E bert sphere initially in 
solid-body rotation ( Bonnor||1956 l. Such analysis allows 
us to compare our results with other numerical simu- 
lations that have considered global collapse and binary 
formation such as |Matsumoto & Hanawal J 2003[) via the 

|T9T)9 



39a I. 



parameters laid out in Tsuribc fc inutsuka] ^ 

We use the collapse calculation of a 10% overdense, 
non - rotati ng Bonnor -Ebert sphere from |Foster fc Cheva-| with 



Fig. 12. — Trajectory of a Bonnor-Ebert sphere through £ — V 
space. The two lines show va lues of /3 = 0.02, 0.08 as defined in 
|Matsumoto fc H anawa ( 2003 1 . Arrows indicate the direction of 
time evolution trom t/tg o = U — 5. tg o is evaluated with respect 
to the central density, and arrows are labelled with the fraction of 
the total Bonnor-Ebert mass which has collapsed up to this point. 
The dotted line shows the fragmentation boundary from Figure [2] 

trajectory of a rotating Bonnor-Ebert sphere through 
£ — r parameter space as a function of the freefall time 
t/tfffi, for two different rotation rates corresponding to 
P = E mt /E gIav = 0.02,0.08. The free fall time is evalu- 
ated with respect to central density. 

The early spike in £ is due to the collapse of the inner 
flattened core at early times. Similarly, the correspond- 
ing decline in V is a result of the mass enclosed increasing 
more rapidly than the inf ailing angular momentum. The 
long period of decreasing £ and constant T arises from 
the balance between larger radii collapsing to contribute 
more angular momentum, and the slow decline of the 
accretion rate. This trajectory may explain several fea- 
tures of the fragmentation seen in|Matsumoto fc Hanawa 



(20031. Although not accounted for in Figure 12 



cores 



lier (1993|, and impose angular momentum on each shell 
to emulate solid body rotation. Figure fl2l shows the 



ligh values of ft have accretion rates supressed at 
early times due to the excess rotational support, while 



those with low (3 collapse at the full rate seen in Foster 
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& Chevalier ( 1993 1. In cores with small (3, the high value 
of £ may drive fragmentation while the disk is young. Al- 
ternatively, for modest values of (3, T may be sufficiently 
low while £ is declining that the disk mass surpasses the 
critical fragmentation threshold, and fragments via the 
so-called satellite formation mechanism. For very large 
values of /3, a core which is only moderately unstable 
will oscillate a nd not collapse as seen in Ma tsumoto fc 
Hanawa| p003| for (3 > 0.3. 



8. DISCUSSION 

We have examined the behavior of gravitationally un- 
stable accretion disks using three-dimensional, AMR nu- 
merical experiments with the code ORION. We charac- 
terize each experiment as a function of two dimensionless 
parameters, £ and T, which arc dimensionless accretion 
rates comparing the infall rate to the disk sound speed 
and orbital period respectively. We find that these two 
global variables can be used to predict disk behavior, 
morphological outcomes, and disk-to-star accretion rates 
and mass ratios. In this first paper in a series we discuss 
the main effects of varying these parameters. Our main 
conclusions are: 

• Disks can process material falling in at up to £ ~ 

2 — 3 without fragmenting. Although increasing T 
stabilizes disks at fixed values of £ those fed at £ > 

3 for many orbits tend to fragment into a multiple 
or binary system. 

• Disks can reach a statistical steady state where 
mass is processed through the disk at a fixed frac- 
tion of the accretion rate onto the disk. The dis- 
crepancy between these two rates, fi, scales with T; 
disks with larger values of T can sustain larger max- 
imum disk masses before becoming unstable. The 
highest disk mass reached in a non-fragmenting sys- 
tem is ju w 0.55 or Af* ~ Md- 

• Gravitational torques can easily produce effective 
accretion rates consistent with a time averaged a ~ 
1. 

• The minimum value of Q at which disks begin to 
fragment is roughly inversely proportional to the 



disk scale height. It is therefore important to con- 
sider not only Q but another dynamical parame- 
ter when predicting fragmentation, at least in disks 
which are not thin and dominated by axisymmetric 
modes. 



• The general disk morphology and multiplicity is 
consistent between isothermal runs and irradiated 
disks with similar effective values of £. 

These conclusions are subject to the qualification that 
fragmentation occurs for lower values of £ as the disk res- 
olution increases, and so it is possible that the location 
of the fragmentation boundary will shift with increasing 
resolution. However we expect that our results are rep- 
resentative of real disks and other numerical simulations 
in so far as they have comparable dynamic range of the 
parameters relevant to fragmentation such as Xj/X. 
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